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Abstract: We describe how regularization of lattice Boltzmann methods can be achieved 
by modifying dissipation. Classes of techniques used to try to improve regularization of 
LBMs include flux limiters, enforcing the exact correct production of entropy and manip- 
ulating non-hydrodynamic modes of the system in relaxation. Each of these techniques 
corresponds to an additional modification of dissipation compared with the standard LBGK 
model. Using some standard ID and 2D benchmarks including the shock tube and lid driven 
cavity, we explore the effectiveness of these classes of methods. 
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1. INTRODUCTION 

Lattice Boltzmann methods (LBM) are a class of discrete computational schemes which can be used to 
simulate hydrodynamics and more|33|. They have been proposed as a discretization of Boltzmann's kinetic 
equation. Instead of fields of moments M, the LBM operates with fields of discrete distributions /. 
All computational methods for continuum dynamics meet some troubles with stability when the gradients of 
^ _ the flows become too sharp. In Computational Fluid Dynamics (CFD) such situations occur when the Mach 

number Ma is not small, or the Reynolds number Re is too large. The possibility for using grid refinement 
r~| \ is bounded by computational time and memory restrictions. Moreover, for nonlinear systems with shocks 

Q-f then grid refinement does not guarantee convergence to the proper solution. Methods of choice to remedy 

this are based on the modification of dissipation with limiters, additional viscosity, and so on |36|. All these 
approaches combine high order methods in relatively quiet regions with low order regularized methods in 
the regions with large gradients. The areas of the high-slope flows are assumed to be small but the loss 
of the order of accuracy in a small region may affect the accuracy in the whole domain because of the 
phenomenon of error propagation. Nevertheless, this loss of accuracy for systems with high gradients seems 
to be unavoidable. 

It is impossible to successfully struggle with some spurious effects without a local decrease in the order 
of accuracy. In a more formal setting, this has been proven. In 1959, Godunov lfT4l proved that a (linear) 
scheme for a partial differintial equation could not, at the same time, be monotone and second-order accurate. 
Hence, we should choose between spurious oscillations in high order non-monotone schemes and additional 
dissipation in first-order schemes. Lax ll23l demonstrated that un-physical dispersive oscillations in areas 
with high slopes are unavoidable due to discretization. For hydrodynamic simulations using the standard 
LBGK model (Sec. „) such oscillations become prevalent especially at high Re and non-small Ma. Levermore 
and Liu used differential approximation to produce the "modulation equation" for the dispersive oscillation 
in the simple initial-value Hopf problem f24l and demonstrated directly how for a nonlinear problem a 
solution of the discretized equation does not converge to the solution of the continuous model with high 
slope when the step /z — > 0. 

Some authors expressed a hope that precisely keeping the entropy balance can make the computation more 
"physical", and that thermodynamics can help to suppress nonphysical effects. Tadmor and Zhong con- 
structed a new family of entropy stable difference schemes which retain the precise entropy decay of the 
Navier-Stokes equations and demonstrated that this precise keeping of the entropy balance does not help to 
avoid the nonphysical dispersive oscillations |.34j. 

To prevent nonphysical oscillations, most upwind schemes employ limiters that reduce the spatial accuracy 
to first order through shock waves. A mixed-order scheme may be defined as a numerical method where 



the formal order of the truncation error varies either spatially, for example, at a shock wave, or for different 
terms in the governing equations, for example, third-order convection with second-order diffusion ||3T1 . 
Several techniques have been proposed to help suppress these pollutive oscillations in LBM, the three which 
we deal with in this work are entropic lattice Boltzman (ELBM), entropic limiters and generalized lattice 
Boltzmann, also known as multiple relaxation time lattice Boltzmann (MRT). Where effective each of these 
techniques corresponds to an additional degree of complexity in the dissipation to the system, above that 
which exists in the LBGK model. 

The Entropic lattice Boltzmann method (ELBM) was invented first in 1998 as a tool for the construction of 
single relaxation time LBM which respect the //-theorem 1 19|. For this purpose, instead of the mirror image 
with a local equilibrium as the reflection center, the entropic involution was proposed, which preserves the 
entropy value. Later, it was called the Karlin-Succi involution iT5\ . In 2000, it was reported that exact 
implementation of the Karlin-Succi involution (which keeps the entropy balance) significantly regularizes 
the post-shock dispersive oscillations |1 1. This regularization seems very surprising, because the entropic 
lattice BGK (ELBGK) model gives a second-order approximation to the Navier-Stokes equation similarly 
to the LBGK model (different proofs of that degree of approximation were given in 1 3 3 1 and 1 4 1 ) . 
Entropic limiters |4| are an example of flux limiter schemes H |2T] |29l , which are invented to combine high 
resolution schemes in areas with smooth fields and first order schemes in areas with sharp gradients. The 
idea of flux limiters can be illustrated by the computation of the flux Fq.i of the conserved quantity u between 
a cell marked by and one of two its neighbour cells marked by ±1: 

Fo,i = (l-0(r))/i7 + 0(r)/^f, (1) 

where /o°*, /o'f^ are low and high resolution scheme fluxes, respectively, r ~ (mq — U-i)/{ui — mq), and 
(^(r) > is a flux limiter function. For r close to 1, the flux limiter function (p(r) should be also close to 1. 
Many flux limiter schemes have been invented during the last two decades |36|. No particular limiter works 
well for all problems, and a choice is usually made on a trial and error basis. Particular examples of the 
limiters we use are introduced in Section |331 

MRT has been developed as a true generalization of the collisions in the lattice Bhatnagar-Gross-Krook 
(LBGK) scheme ifTTl l22l from a one parameter diagonal relaxation matrix, to a general linear operation 
with more free parameters, the number of which is dependent on the particular discrete velocity set used and 
the number of conserved macroscopic variables. Different variants of MRT have been shown to improve 
accuracy and stability, including in our benchmark examples ll26l in comparison with the standard LBGK 
systems. 

The lattice Boltzmann paradigm is now mature, and explanations for some of its successes are available. 
However, in its applications it approaches the boundaries of the applicability and need special additional 
tools to extend the area of applications. It is well-understood that near to shocks, for instance, special and 
specific attention must be paid to avoid unphysical effects. In this paper then, we will discuss a variety of 
add-ons for LBM and apply them to a variety of standard ID and 2D problems to test their effectiveness. In 
particular we will describe a family of entropic filters and show that we can use them to signficantly expand 
the effective range of operation of the LBM. 

2, BACKGROUND 

Lattice Boltzmann methods can be derived independently by a discretization Boltzmann's equation for ki- 
netic transport or by naively creating a discrete scheme which matches moments with the Maxwellian dis- 
tribution up to some finite order. 

In each case the final discrete algorithm consists of two alternating steps, advection and collision, which 
are applied to m single particle distribution functions fi = fi{x,t), (/=!.. .m), each of which corresponds 
with a discrete velocity vector v,-, (/ = 1 . ..m). The values /,■ are also sometimes known as populations or 
densities as they can be thought of as representative of the densities of particles moving in the direction of 
the corresponding discrete velocities. 

The advection operation is simply free flight for the discrete time step e in the direction of the corresponding 
velocity vector. 



fi{x,t + e)=fi{x-evi,t). (2) 

The collision operation is instantaneous and can be different for each distribution function but depends on 
every distribution function, this might be written, 

Mx)^F,{{Mx)}). (3) 

In order to have a slightly more compact notation we can write these operations in vector form, in the below 
equation it should be inferred that the /th distribution function is advecting along its corresponding discrete 
velocity, 

f(x,f + e) = f(x-ev,-,f), (4) 

f(x)^F(f(x)). (5) 

To transform our vector of microscopic variables at a point in space f(x) to a vector of macroscopic vari- 
ables M(x) we use a vector of linear functions M(x) = m(f(x)). In the athermal hydrodynamic systems we 
consider in this work the momemts are density p and momentum density pu, {p,pu}(x) = Y.i{^i'^i}.fi{'^)- 
These macroscopic moments are conserved by the collision operation, m(f) = m{F{f)). 
The simplest and most common choice for the collision operation F is the Bhatnagar-Gross-Krook(BGK) 17] 
|9] [l7| i33J operator with over-relaxation 

F(f)=f+aj3(n-f). (6) 

For the standard LBGK method a~2 and j3 G [0, 1] (usually, j3 G [1/2, 1]) is the over-relaxation coefficient 
used to control viscosity. For j3 = 1 /2 the collision operator returns the local equilibrium F'i and j3 = 1 (the 
mirror reflection) returns the collision for a liquid at the zero viscosity limit. The definition of F'^ defines the 
dynamics of the system, often it chosen as an approximation to the continuous Maxwellian distribution. An 
equilibrium can also be independently derived by constructing a discrete system which matches moments 
of the Maxwellian up to some finite order. For hydrodynamic systems often this finite order is chosen to 
be 2, as this is sufficient to accurately replicate the Euler(non dissipative) component of the Navier Stokes 
equations. For a dissipative fluid with viscosity v the parameter )3 is chosen by )3 = e/(2v + e). 
Each of the techniques we test in this paper can be introduced as developments of the generic LBGK system 
and such a presentation follows in the next sections. 

3. ENTROPIC LATTICE BOLTZMANN 



3.1. LBM with H theorem 

In the continuous case the Maxwellian distribution maximizes entropy, as measured by the Boltzmann H 
function, and therefore also has zero entropy production. In the context of lattice Boltzmann methods a 
discrete form of the //-theorem has been suggested as a way to introduce thermodynamic control to the 
system ITSlfTi. 

From this perspective the goal is to find an equilibrium state equivalent to the Maxwellian in the continuum 
which will similarly maximize entropy. Before the equilibrium can be found an appropriate H function must 
be known for a given lattice. These functions have been constructed in a lattice dependent fashion in IITSjl . 
and H = —S with S from Q is an example of a // function constructed in this way. 

One way to implement an ELBM is as a variation on the LBGK, known as the ELBGK In this case a 
is varied to ensure a constant entropy condition according to the discrete //-theorem. In general the entropy 
function is based upon the lattice and cannot always be found explicitly. However for some examples 
such as the simple one dimensional lattice with velocities v = (— c,0,c) and corresponding populations 
f = (/-,/o,/+) an explicit Boltzmann style entropy function is known IfTSl : 



5(f) = -/_log(/_) -/olog(/o/4) -/+log(/+). 



(7) 



With knowledge of such a function a is found as the non-trivial root of the equation 



5(f) = 5(f+a(r-f)). 



(8) 



The trivial root a = returns the entropy value of the original populations. ELBGK then finds the non-trivial 
a such that (jSj holds. This version of the BGK collision one calls entropic BGK (or EBGK) collision. A 
solution of dHJ must be found at every time step and lattice site. Entropic equilibria (also derived from the 
//-theorem) are always used for ELBGK. 

3.2. ELBM algorithm and additional dissipation 

The definition of ELBM for a given entropy equation ([8]) is incomplete. First of all, it is possible that the 
non-trivial solution does not exist. Moreover, for most of the known entropies (like the perfect entropy [18 ]) 
there always exist such / that the equation (|8]l for the ELBM collision has no non-trivial solutions. These / 
should be sufficiently far from equilibrium. For completeness, every user of ELBM should define collisions 
when the non-trivial root of dH) does not exist. We know and tried two rules for this situation: 

1. The most radical approach gives the the Ehrenfest rule lfT6l l3l: "if the solution does not exist then go 
to equilibrium", i.e. if the solution does not exists then take a = 1. 

2. The most gentle solution gives the "positivity rule" |I3] |25] [35] [32l : to take the maximal value of a 
that guarantees /, + a{f * — /,) > for all /. 

In general, the Ehrenfest rule prescribes to send the most non-equilibrium sites to equilibrium and the posi- 
tivity rule is applied for any LBM as a recommendation to substitute the non-positive vectors f by the closest 
non-negative on the interval of the straight line [f, T] that connects f to equilibrium. These rules give the 
examples of the pointwise LBM limiter and we discuss them separately. 

By its nature, the ELBM adds more dissipation than the positivity rule when the non-trivial root of (O does 
not exist. It does not always keep the entropy balance but increases dissipation for highly nonequilibrium 
sites. 

3.3. Numerical method for solution of the ELBM equation 

Another source for additional dissipation in the ELBM may be the numerical method used for the solution of 
([8]). For the full description of ELBM we have to select a numerical method for this equation. This method 
has to have an uniform accuracy in the wide range of parameters, for all possible deviation from equilibrium 
(distribution of these deviations has "heavy tails" |5|). 

In order to investigate the stabilization properties of ELBGK it is necessary to craft a numerical method 
capable of finding the non-trivial root in (O. In this section we fix the population vectors f and f*, and are 
concerned only with this root finding algorithm. We recast (O as a function of a only; 



In this setting we attempt to find the non-trivial root r of (|9ll such that Sf{r) — 0. It should be noted that as we 
search for r numerically we should always take care that the approximation we use is less than r itself. An 
upper approximation could result in negative entropy production. A simple algorithm for finding the roots 
of a concave function, based on local quadratic approximations to the target function, has cubic convergence 
order. Assume that we are operating in a neighbourhood r G A^, in which S'^ is negative (as well of course S" 
is negative). At each iteration the new estimate for r is the greater root of the parabola P, the second order 
Taylor polynomial at the current estimate. Analogously to the case for Newton iteration, the constant in the 
estimate is the ratio of third and first derivatives in the interval of iteration: 



Sf{a)^S{{+a{n-{))-S{{). 



(9) 



\{r-a„+i)\ < C\a„-r\^, 




where a„ is the evaluation of r on the nth iteration. 

We use a Newton step to estimate the accuracy of the method at each iteration: because of the concavity of 
S 

\cc„-r\<\Sf{a„)/S'f{a„)\. (10) 

In fact we use a convergence criteria based not solely on a but on a||f* — f||, this has the intuitive appeal 
that in the case where the populations are close to the local equilibrium AS = S{f*) — 5(f) will be small and 
a very precise estimate of a is unnecessary. We have some freedom in the choice of the norm used and we 
select between the standard Li norm and the entropic norm. The entropic norm is defined as 

||fq_f||p,^_((fq_f),Z)25|p,(n-f)), 

where is the second differential of entropy at point F'', and {x,y) is the standard scalar product. 

The final root finding algorithm then is beginning with the LBGK estimate ;ico = 2 to iterate using the roots 
of successive parabolas. We stop the method at the point, 

|a„-rH|r-f|| <e. (11) 

To ensure that we use an estimate that is less than the root, at the point where the method has converged we 
check the sign of Sf{a„). If Sf{an) > then we have achieved a lower estimate, if Sf{a„) < we correct 
the estimate to the other side of the root with a double length Newton step, 

Sf(a„) 

^"-^"-'^y '''' 

At each time step before we begin root finding we eliminate all sites with A5 < 10^'^. For these sites we 
make a simple LBGK step. At such sites we find that round off error in the calculation of Sf by solution of 
equation (|8]l can result in the root of the parabola becoming imaginary. In such cases a mirror image given 
by LBGK is effectively indistinct from the exact ELBGK collision. In the numerical examples given in this 
work the case where the non-trivial root of the entropy parabola does not exist was not encountered. 

4. ENTROPIC FILTERING 

All the specific LBM limiters |I5]|29] are based on a representation of distributions / in the form: 

f=feq+||f_feq||^5^, (13) 

where f is the corresponding quasiequilibrium (conditional equilibrium) for given moments M, f — F'^ is the 
nonequilibrium "part" of the distribution, which is represented in the form "norm x direction" and ||/ — /*|| 
is the norm of that nonequilibrium component (usually this is the entropic norm). 

All limiters we use change the norm of the nonequilibrium component f — pi, but do not touch its direction 
or the quasiequilibrium. In particular, limiters do not change the macroscopic variables, because moments 
for f and 1*'^ coincide. These limiters are transformations of the form 

f^->f'^ + X (f-fi) (14) 

with (/) > 0. If f — F'' is too big, then the limiter should decrease its norm. 

For the first example of the realization of this pointwise filtering we use the kinetic idea of the positivity rule, 
the prescription is simple |I3] |25] |35] |32 1 : to substitute nonpositive F (f ) by the closest nonnegative state that 
belongs to the straight line 

|Af+(l -A)fi|A g]r| (15) 

defined by the two points, f and the corresponding quasiequilibrium [(T)| This operation is to be applied 
pointwise, at points of the lattice where positivity is violated. This technique preserves the positivity of 
populations, but can affect the accuracy of the approximation. This rule is necessary for ELBM when the 
positive "mirror state" with the same entropy as f does not exists on the straight line (flST l. 




Fig. (1): Positivity rule in action. The motions stops at the positivity boundary. 



The positivity rule measures the deviation / — /* by a binary measure: if all components of this vector 
f — f* are non-negative then it is not too large. If some of them are negative that this deviation is too large 
and needs corrections. 

To construct pointwise flux limiters for LBM, based on dissipation, the entropic approach remains very 
convenient. The local nonequilibrium entropy for each site is defined 

A5(f) :=5(fi)-5(f). (16) 

The positivity limiter was targeted wherever population densities became negative. Entropic limiters are 
targeted wherever non-equilibrium entropy becomes large. 

The first limiter is Ehrenfests' regularisation ||4]|3], it provides "entropy trimming": we monitor local de- 
viation of f from the corresponding quasiequilibrium, and when A5(f)(x) exceeds a pre-specified threshold 
value 5, perform local Ehrenfests' steps to the corresponding equilibrium: f H> at those points. 
Not all lattice Boltzmann models are entropic, and an important question arises: "how can we create 
nonequilibrium entropy limiters for LBM with non-entropic (quasi)equilibria?". We propose a solution 
of this problem based on the discrete Kullback entropy [20.1 : 

5K(f) = -E/,ln(^-^). (17) 

For entropic quasiequilibria with perfect entropy the discrete Kullback entropy gives the same A5: — 5k(/) = 
A5(/). Let the discrete entropy have the standard form for an ideal (perfect) mixture l,18j : 

In quadratic approximation, 

-5K(f)-E/'ln(;^) «E^^^^#^- (18) 
If we define f as the conditional entropy maximum for given Mj = Y.k^'^ikfk^ then 

J 

where pij{M) are the Lagrange multipliers (or "potentials"). For this entropy and conditional equilibrium 
we find 

A5 = 5(fi)-5(f) = E/;in(^A^ =_5(,(f), (19) 
if f and F'^ have the same moments, m(f) = m(F'i). 

In what follows, A5 is the Kullback distance — 5K(f^'') (O for general (positive) quasiequilibria 1*1, or 
simply 5(1*'') — 5(f) for entropic quasiequilibria (or second approximations for these quantities (fTsTl). 



So that the Ehrenfests' steps are not allowed to degrade the accuracy of LBGK it is pertinent to select the k 
sites with highest AS > d. An a posteriori estimate of added dissipation could easily be performed by the 
analysis of entropy production in the Ehrenfests' steps. Numerical experiments show (see, e.g., |4, 3|) that 
even a small number of such steps drastically improves stability. 

The positivity rule and Ehrenfests' regularisation provide rare, intense and localised corrections. Of course, 
it is easy and also computationally cheap to organise more gentle transformations with a smooth shift 
of highly nonequilibrium states to quasiequilibrium. The following regularisation transformation with a 
smooths function distributes its action smoothly: 

fH^fi + 0(A5(f))(f-n). (20) 

The choice of function (j) is highly ambiguous, for example, = 1/(1 + aAS'') for some a > and A: > 0. 
There are two significantly different choices: (i) ensemble-independent (j) (i.e., the value of (j) depends on 
local value of AS only) and (ii) ensemble-dependent (j), for example 

^^^'"'^ l + {AS/{aE{AS))r ' ^^^^ 

where E(A5) is the average value of AS in the computational area, k>l, and o: > 1 . For small AS, (p (AS) ~ 1 
and for AS 3> aE(A5) it tends to ^/aE{AS)/AS. It is easy to select an ensemble-dependent (p with control 
of total additional dissipation. 



4.1. Monotonic and double monotonic limiters 

Two monotonicity properties are important in the theory of nonequilibrium entropy limiters: 

1. A limiter should move the distribution to equilibrium: in all cases of ( fT4] i < (j) < I. This is the 
dissipativity condition which means that limiters never produce negative entropy. 

2. A limiter should not change the order of states on the line: if for two distributions with the same 
moments, f and f , f — Fi = x(f — f) and A5'(f) > A5'(f ) before the limiter transformation, then the 
same inequality should hold after the limiter transformation too. For example, for the limiter ( l20b it 
means that A5(F1 + x0(A5(pi +x(f — f"i))(f — fi)) is a monotonically increasing function of x > 0. 

In quadratic approximation, 

A5(f'J+x(f-fi))=x2A5(f), 
A5 (n + ( A5 (f + X (f - n ) (f - f'' ) ) = x2 2 (^2^5 (f ) ) , 

and the second monotonicity condition transforms into the following requirement: y<p {y^s) is a monotoni- 
cally increasing (not decreasing) function of y > for any 5 > 0. 

If a limiter satisfies both monotonicity conditions, we call it "double monotonic". For example, Ehrenfests' 
regularisation satisfies the first monotonicity condition, but violates the second one. The limiter ( |2TI ) violates 
the first condition for small AS, but is dissipative and satisfies the second one in quadratic approximation 
for large AS. The limiter with = 1/(1 + aAS'^) always satisfies the first monotonicity condition, violates 
the second if A: > 1 /2, and is double monotonic (in quadratic approximation for the second condition), if 
< A; < 1/2. The threshold limiter ( |26] | is also double monotonic. 

For smooth functions, the condition of double monotonicity (in quadratic approximation) is equivalent to 
the system of differential inequalities: 

<p{x)+2x<p'{x) >0; 

0'(jc) <0. 

The initial condition (0) = 1 means that in the limit A5 ^ limiters do not affect the flow. Following these 
inequalities we can write: 2x0 '(x) = — Tj(x)0(x), where < Tj(x) < 1. The solution of these inequalities 
with initial condition is 

0W=exp(^-^^'^d;t), (22) 



if the integral on the right-hand side exists. This is a general solution for double monotonic limiters (in the 
second approximation for entropy). If rj (x) is the Heaviside step function, Tj (x) — H{x — ASi) with threshold 
value ASt, then the general solution (l22l l gives us the threshold limiter. If, for example, ri{x) =x*^/(A5f 
then 

_ j_ 

/ a- \ Ik 

+ . (23) 

This special form of limiter function is attractive because for small x it gives 

1 x'' 
2fcA5f 



Thus, the limiter does not affect the motion up to the (^-|- l)st order, and the macroscopic equations coincide 
with the macroscopic equations for LBM without limiters up to the {k+ l)st order in powers of deviation 
from quasiequilibrium. Furthermore, for large x we get the ^th order approximation to the threshold limiter 




Of course, it is not forbidden to use any type of limiters under the local and global control of dissipation, but 
double monotonic limiters provide some natural properties automatically, without additional care. 

4.2. Monitoring total dissipation 

For given j3, the entropy production in one LBGK step in quadratic approximation for AS is: 

5LBGK5«[l-(2j3-l)2]^A5(x), 

X 

where x is the grid point, A5'(x) is nonequilibrium entropy ( fT6b at point x, 5lbgk5' is the total entropy 
production in a single LBGK step. It would be desirable if the total entropy production for the limiter 5\imS 
was small relative to 5lbgk5': 

diimS < 5o5lbgk5. (24) 

A simple ensemble-dependent limiter (perhaps, the simplest one) for a given 5o operates as follows. Let 
us collect the histogram of the A5(x) distribution, and estimate the distribution density, p{AS). We have to 
estimate a value ASq that satisfies the following equation: 

/*oo 

p(A5)(A5-A5o)dA5 = 5o[l-(2j3-l)2] / p{AS)ASdAS. (25) 

ASo ^0 

In order not to affect distributions with a small expectation of AS, we choose a threshold A5t = max{A5o, 5}, 
where 5 is some predefined value (as in the Ehrenfests' regularization). For states at sites with AS > ASt 
we provide homothety with equilibrium center and coefficient y/ ASt /AS (in quadratic approximation for 
nonequilibrium entropy): 

f(x) ^f=''(x) + y^(f(x) -n(x)). (26) 

To avoid the change of accuracy order "on average", the number of sites with this step should be < {Nh /L) 
where is the total number of sites, h is the step of the space discretization and L is the macroscopic charac- 
teristic length. But this rough estimate of accuracy across the system might be destroyed by a concentration 
of Ehrenfests' steps in the most nonequilibrium areas, for example, in boundary layers. In that case, instead 
of the total number of sites in €>'{Nh/L) we should take the number of sites in a specific region. The 
effects of such concentration could be analysed a posteriori. 



4.3. Median entropy filter 



The Ehrenfest step described above provides pointwise correction of nonequilibrium entropy at the "most 
nonequilibrium" points. Due to the pointwise nature, the technique does not introduce any nonisotropic 
effects, and provides some other benefits. But if we involve local structure, we can correct local non- 
monotone irregularities without touching regular fragments. For example, we can discuss monotone in- 
crease or decrease of nonequilibrium entropy as regular fragments and concentrate our efforts on reduction 
of "speckle noise" or "salt and pepper noise". This approach allows us to use the accessible resource of 
entropy change ( |24| ) more thriftily. Salt and pepper noise is a form of noise typically observed in images. 
It represents itself as randomly occurring white (maximal brightness) and black pixels. For this kind of 
noise, conventional low-pass filtering, e.g., mean filtering or Gaussian smoothing is unsuccessful because 
the perturbed pixel value can vary significantly both from the original and mean value. For this type of noise, 
median filtering is a common and effective noise reduction method. Median filtering is a common step in 
image processing Il30ll for the smoothing of signals and the suppression of impulse noise with preservation 
of edges. 

The median is a more robust average than the mean (or the weighted mean) and so a single very unrepresen- 
tative value in a neighbourhood will not affect the median value significantly. Hence, we suppose that the 
median entropy filter will work better than entropy convolution filters. 

For the nonequilibrium entropy field, the median filter considers each site in turn and looks at its nearby 
neighbours. It replaces the nonequilibrium entropy value AS at the point with the median of those values 
AS'niedi then updates / by the transformation ( |26] ) with the homothety coefficient y/ ASi-ned/^S. The median, 
A5nied, is calculated by first sorting all the values from the surrounding neighbourhood into numerical order 
and then replacing that being considered with the middle value. For example, if a point has 3 nearest 
neighbours including itself, then after sorting we have 3 values AS: ASi < AS2 < AS3. The median value is 
^•Smed = ^'^2- For 9 nearest neighbours (including itself) we have after sorting A^med = ^S^. For 27 nearest 
neighbours A5med = A5i4. 

We accept only dissipative corrections (those resulting in a decrease of AS, A5med < ^S) because of the 
second law of thermodynamics. The analogue of (1251 ) is also useful for the acceptance of the most significant 
corrections. In "salt and pepper" terms, we correct the salt (where AS exceeds the median value) and do not 
touch the pepper. 

4.4. General nonlocal filters 

The separation of f in equilibrium and nonequilibrium parts (fTJb allows one to use any nonlocal filtering 
procedure. Let — f — F^. The values of moments for f and f"* coincide, hence we can apply any 
transformation of the form 

fneq(x)^£rf(y)fneq(x + y) 

_v 

for any family of vectors y that shift the grid into itself and any coefficients d{y). If we apply this trans- 
formation, the macroscopic variables do not change but their time derivatives may change. We can control 
the values of some higher moments in order not to perturb significantly some macroscopic parameters, the 
shear viscosity, for example LIOJ . Several local (but not pointwise) filters of this type have been proposed 
and tested recently ll29l . 

5. MULTIPLE RELAXATION TIMES 

The MRT lattice Boltzmann system Il22l l26l [TTII generalizes the BGK collision into a more general linear 
transformation of the population functions, 

F(f) =f + A(fi-f), (27) 

where A is a square matrix of size m. The use of this more general operator allows more different parameters 
to be used within the collision, to manipulate different physical properties, or for stability purposes. 
To facilitate this a change of basis matrix can be used to switch the space of the collision to the moment 
space. Since the moment space of the system may be several dimensions smaller than the population space. 



to complete the basis linear combinations of higher order polynomials of the discrete velocity vectors may 
be used. For our later experiments we will use the D2Q9 system, we should select a particular enumeration 
of the discrete velocity vectors for the system, the zero velocity is numbered one and the positive x velocity 
is numbered 2, the remainder are numbered clockwise from this system, 

7 8 9 

6 12. (28) 
5 4 3 

The change of basis matrix given in some of the literature ll22l is chosen to represent specific macroscopic 
quantities and in our velocity system is as follows. 



/ 


1 


1 


1 


1 


1 


1 


1 


1 


1 \ 




-4 


-1 


2 


-1 


2 


-1 


2 


-1 


2 




4 


-2 


1 


-2 


1 


-2 


1 


-2 


1 







1 


1 





-1 


-1 


-1 





1 







-2 


1 





-1 


2 


-1 





1 










-1 


-1 


-1 





1 


1 


1 










-1 


2 


-1 





1 


-2 


1 







1 





-1 





1 





-1 





V 








-1 





1 





-1 





1 / 



(29) 



As an alternative we could complete the basis simply using higher powers of the velocity vectors, the basis 
would be 1, v'l-, V,., + v^, — v^, VjcV-y, vjvy, V;cVy, v^Vy, in our velocity system then this change of basis matrix 
is as follows. 
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When utilized properly any basis should be equivalent (although with different rates). In any case in this 
system there are 3 conserved moments and 9 population functions. Altogether then there are 6 degrees of 
freedom in relaxation in this system. We need some of these degrees of freedom to implement hydrodynamic 
rates such as shear viscosity or to force isotropy and some are 'spare' and can be manipulated to improve 
accuracy or stability. Typically these spare relaxation modes are sent closer to equilibrium than the standard 
BGK relaxation rate. This is in effect an additional contraction in the finite dimensional non-equilibrium 
population function space, corresponding to an increase in dissipation. 

Once in moment space we can apply a diagonal relaxation matrix Ci to the populations and then the inverse 
moment transformation matrix Mj"' to switch back into population space, altogether A i = M^^CiMi. If 
we use the standard athermal polynomial equilibria then three entries on the diagonal of Ci are actually not 
important as the moments will be automatically conserved since m(f) = m(l*i), for simplicity we set them 
equal to or 1 to reduce the complexity of the terms in the collision matrix. There are 6 more parameters 
on the diagonal matrix C which we can set. Three of these correspond to second order moments, one each 
is required for shear and bulk viscosity which are called Se and Sy respectively and one for isotropy. Two 
correspond to third order moments, one gives a relaxation rate and again one is needed for isotropy. 
Finally one is used to give a relaxation rate for the single fourth order moment. We have then in total four 



relaxation parameters which appear on the diagonal matrix in the following form: 
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(31) 



Apart from the parameter Sy which is used to control shear viscosity, in an incompressible system the other 
properties can be varied to improve accuracy and stability. In particular, there exists a variant of MRT known 
as TRT (two relaxation time) lfT3l where the relaxation rates ie,i£ are made equal to Sy. In a system with 
boundaries the final rate is calculated Sg — ^{2 — Sy) / — Sy), this is done, in particular, to combat numerical 
slip on the boundaries of the system. 

We should say that in some of the literature regarding MRT the equilibrium is actually built in moment 
space, that is the collision operation would be written. 



F(f) f + M-^Ci(m'=i -Mif). 



(32) 



This could be done to increase efficiency, depending on the implementation, however each moment equi- 
librium m'^i has an equivalent population space equilibrium f^'' — M^^af^, so the results of implementing 
either system should be the same up to rounding error. 

We can also conceive of using an MRT type collision as a limiter, that is to apply it only on a small number 
of points on the lattice where non equilibrium entropy passes a certain threshold. This answers a criticism 
of the single relaxation time limiters that they fail to preserve dissipation on physical modes. As well as 
using the standard MRT form given above we can build an MRT limiter using the alternative change of basis 
matrix M2. 

The limiter in this case is based on the idea of sending every mode except shear viscosity directly to equi- 
librium again the complete relaxation matrix is given by A2 = M^^C2M2 where. 
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(33) 



This could be considered a very aggressive form of the MRT which maximizes regularization on every mode 
except shear viscosity and would not be appropriate for general use in a system, especially as most systems 
violate the incompressibility assumption and hence bulk viscosity is not small. The advantage of using the 
different change of basis matrix is that the complete collision matrix A2 is relatively sparse with just twelve 
off diagonal elements and hence is easy to implement and not too expensive to compute with. 

6. ID SHOCK TUBE 

A standard experiment for the testing of LBMs is the one-dimensional shock tube problem. The lattice 
velocities used are v = (—1,0, 1), so that space shifts of the velocities give lattice sites separated by the unit 
distance. 800 lattice sites are used and are initialized with the density distribution 



p{x) = 



1, 

0.5, 



l<x< 400, 
401<jc < 800. 



Initially all velocities are set to zero. We compare the ELBGK equipped with the parabola based root 
finding algorithm using the entropic norm with the standard LBGK method using both standard polynomial 
and entropic equilibria. The polynomial equilibria are given in I7l l33i : 

f* P (-1 1 , ^ 2\ r* 2p / 3m2\ 

/--■^(i-3m+3m ), fi) = — y-—y 

/+ = |(l+3« + 3«'). 

The entropic equilibria also used by the ELBGK are available explicitly as the maximum of the entropy 
function O, 

/* =£(_3„_l+2v/l+3«2), /* = ^(2-v/1+3m2), 
/; = |(3m-1+2\/1+3m2). 

Now following the prescription fromm Sec. □ the governing equations for the simulation are 

/_(x,f+l)=/_(x+l,f)+aj3(/:(x+l,f)-/-(x+l,0), 
/o (x, f + 1 ) =/o (x, f ) + aj3 (/o* (x, f ) - /o (x, f )) , 

/+(x,f+l)=/+(x-l,f)+ai8(/;(x-l,f)-/+(x-l,f)). 
From this experiment we observe no benefit in terms of regularization in using the ELBGK rather than the 
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Fig. (2): Density profile of the simulation of the shock tube problem following 400 time steps using (a) 
LBGK with polynomial equilibria [v = (1/3) • 10"']; (b) LBGK with entropic equilibria [v = (1/3) • 10"']; 
(c) ELBGK [V = (1 /3) • 10"']; (d) LBGK with polynomial equilibria [v = lO"*^]; (e) LBGK with entropic 
equilibria [v = lO"*^]; (f) ELBGK [v = lO-**]. 



standard LBGK method (Fig. |(2)| l. In both the medium and low viscosity regimes ELBGK does not supress 
the spurious oscillations found in the standard LBGK method. The observation is in full agreement with the 
Tadmor and Zhong Il34l experiments for schemes with precise entropy balance. 



Entropy balance gives a nice additional possibility to monitor the accuracy and the basic physics but does 
not give an omnipotent tool for regularization. 

7. 2D SHEAR DECAY 

For the second test we use a simple test proposed to measure the observable viscosity of a lattice Boltz- 
mann simulation to validate the shear viscosity production of the MRT models. We take the 2D isother- 
mal nine-velocity model with standard polynomial equilibria. Our computational domain will a square 
which we discretize with L + I x L + I uniformly spaced points and periodic boundary conditions. The 
initial condition is p{x,y) = l,Ux{x,y) = and Uy{x,y) — uosin{2nx/L), with mq = 0.05. The exact ve- 
locity solution to this problem is an exponential decay of the initial condition: Ux{x,y,t) = Q,Uy{x,y,t) = 
Moexp(— AMo?/ReL) sin(2;rjic/L), where A is some constant and Re = uqL/v is the Reynolds number of the 
flow. Here, v is the theoretical shear viscosity of the fluid due to the relaxation parameters of the collision 
operation. 

Now, we simulate the flow over L/uq time steps and measure the constant A from the numerical solution. 
We do this for LBGK, the 'aggressive' MRT with collision matrix A2 and the MRT system with collision 
matrix A 1 with additional parameters Sg = 1.64,i£ = 1.54,i^ = 3(2 — iv)/(3 — Sy) Il22ll . The shear viscosity 
relaxation parameter Sv is varied to give different viscosities and therefore Reynolds numbers for L = 50 and 
for L= 100. 




Fig. (3): Observed A in the shear decay experiment 

From Figure [(3)1 it can be seen that in all the systems, for increasing theoretical Reynolds numbers (decreas- 
ing viscosity coefficient) following a certain point numerical dissipation due to the lattice begins growing. 
The most important observation from this system is that our MRT models do indeed produce shear viscosity 
at almost the same rate as the BGK model. 

In particular, the bulk viscosity in this system is zero, so all dissipation is given by shear viscosity or higher 
order modes. Since the space derivatives of the velocity modes are well bounded, as is the magnitude of the 
velocity itself, the proper asymptotic decay of higher order modes is observed and the varying higher order 
relaxation coefficients have only a very marginal effect. 

We should reiterate that we while we use the 'aggressive' MRT across the system, this is only appropriate 
as bulk viscosity is zero and in fact the selection of the bulk viscosity coefficient makes no difference in this 
example. This example is a special case in this regard. 



8. LID DRIVEN CAVITY 



8.1. Stability 

Our next 2D example is the benchmark 2D lid driven cavity. In this case this is a square system of side 
length 129. Bounce back boundary conditions are used and the top boundary imposes a constant velocity of 
Mwaii =0.1. For a variety of Reynolds numbers we run experiments for up to 10000000 time steps and check 
which methods have remained stable up until that time step. 

The methods which we test are the standard BGK system, the BGK system equipped with Ehrenfest steps, 
the BGK system equipped with the MRT limiter, the TRT system, an MRT system with the TRT relaxation 
rate for the third order moment and the other rates Se = l.64,S£ = 1.54 and finally an MRT system which 
we call MRTl with the rates = l.9se = 1.64, Sg = 1.54 |26|. In each case of the system equipped with 
limiters the maximum number of sites where the limiter is used is 9. 

All methods are equipped with the standard 2nd order compressible quasi-equilibrium, which is available as 
the product of the ID equilibria l4.4.l 

When calculating the stream functions of the final states of these simulations we use Simpson integration in 
first the X and then y directions. 

Additionally we measure Enstrophy S" in each system over time. Enstrophy is calculated as the sum of 
vorticity squared across the system, normalized by the number of lattice sites. This statistic is useful as 
vorticity is theoretically only dissipated due to shear viscosity, at the same time in the lid driven system 
vorticity is produced by the moving boundary. For these systems becomes constant as the vorticity field 
becomes steady. The value of this constant indicates where the balance between dissipation and production 
of vorticity is found. The lower the final value of S' the more dissipation produced in the system. 
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Fig. (4): Contour plots of stream functions of A: BGK, B: BGK + Ehrenfest Steps, C: BGK + MRT Limiter, 
D: TRT, E: MRT, F: MRTl following 10000000 time steps at RelOOO. 

All of these systems are stable for RelOOO, the contour plots of the final state are given in Figure [(4)| and 
there appears only small differences. We calculate the average enstrophy in each system and plot it as a 
function of time in Figure [(5)| We can see that in the different systems that enstrophy and hence dissipation 
varies. Compared with the BGK system all the other systems except MRTl exhibit a lower level of enstrophy 
indicating a higher rate of dissipation. For MRTl the fixed relaxation rate of the third order mode is actually 
less dissipative than the BGK relaxation rate for this Reynolds number, hence the increased enstrophy. An 
artifact of using the pointwise filtering techniques is that they introduce small scale local oscillations in 
the modes that they regularize, therefore the system seems not to be asymptotically stable. This might be 
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Fig. (5): Enstrophy in the RelOOO systems during the final 2 • 10 time steps 



remedied by increasing the threshold of AS below which no regularization is performed. Nevertheless in 
these experiments after sufficient time the enstrophy values remain within a small enough boundary for the 
results to be useful. 




Fig. (6): Contour plots of stream functions of A: BGK + Ehrenfest Steps, B: BGK + MRT Limiter, C: TRT, 
D: MRT 1, following 10000000 time steps at Re2500. 
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Fig. (7): Enstrophy in the Re2500 systems during the final 2 • 10^ time steps 



The next Reynolds number we choose is Re2500. Only 4 of the original 6 systems complete the full number 
of time steps for this Reynolds number, the contour plots of the final stream functions are given in Figure 
1(6)1 Of the systems which did not complete the simulation it should be said that the MRT system survived 
a few 10s of thousands of time steps while the BGK system diverged almost immediately, indicating that it 
does provide stability benefits which are not apparent at the coarse granularity of Reynolds numbers used in 
this study. One feature to observe in the stream function plots is the absence of an upper left vortex in the 
Ehrenfest limiter. This system selects the "most non-equilibrium" sites to apply the filter. These typically 
occur near the corners of the moving lid. It seems here that the local increase in shear viscosity is enough 
to prevent this vortex forming. This problem does not seem to affect the MRT limiter which preserves the 
correct production of shear viscosity. 

Again we check the enstrophy of the systems and give the results for the final timesteps in Figure [(7)1 Due 
to the failure of the BGK system to complete this simulation there is no "standard" result to compare the 
improved methods with. The surviving methods maintain their relative positions with respect to enstrophy 
production. 
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Fig. (8): Contour plots of stream functions of A: BGK + Ehrenfest Steps, B: MRT 1, following 10000000 
time steps at ReSOOO. 

For the theoretical Reynolds number of 5000 only two systems remain, their streamfunction plots are given 
in Figure [(8)l At this Reynolds number the upper left vortex has appeared in the Ehrenfest limited simulation, 
however a new discrepancy has arisen. The lower right corner exhibits a very low level of streaming. 
In Figure [(9)] the enstrophy during the final parts of the simulation is given. We note that for the first time the 



- B - BGK+Ehr 
MRT1 



Fig. (9): Enstrophy in the ReSOOO systems during the final 2 • 10^ time steps 




Fig. (10): Contour plots of stream functions of A: Re7500 and B: RelOOOO BGK + Ehrenfest sys- 
tems,following 10000000 time steps at ReSOOO. 

For the final two Reynolds numbers we use, 7500 and 10000, only the BGK system with the Ehrenfest 
limiter completes the simulation. The corresponding streamfunction plots are given in Figure [(T0)| and they 
exhibit multiple vortices in the corners of the domain. 

The enstrophy plots are given in Figure [(TT)] as the theoretical Reynolds number increases so does the level 
of enstrophy. 

8.2. First Hopf Bifurcation 

As Reynolds number increases the flow in the cavity is no longer steady and a more complicated flow pattern 
emerges. On the way to a fully developed turbulent flow, the lid-driven cavity flow is known to undergo a 
series of period doubling Hopf bifurcations. 

A survey of available literature reveals that the precise value of Re at which the first Hopf bifurcation occurs 
is somewhat contentious, with most current studies (all of which are for incompressible flow) ranging from 
around Re = 7400-8500 161 l27l l28l . Here, we do not intend to give a precise value because it is a well 
observed grid effect that the critical Reynolds number increases (shifts to the right) with refinement (see. 
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Fig. (11): Enstrophy in the Re7500 and RelOOOO BGK + Ehrenfest systems during the final 2 10 time steps 



e.g.. Fig. 3 in [28 1). Rather, we will be content to localise the first bifurcation and, in doing so, demonstrate 
that limiters are capable of regularising without effecting fundamental flow features. 

To localise the first bifurcation we take the following algorithmic approach. Entropic equilibria are in use. 
The initial uniform fluid density profile is p = 1.0 and the velocity of the lid is uq = 1/10 (in lattice units). 
We record the unsteady velocity data at a single control point with coordinates (L/16, 13L/16) and run 
the simulation for 5000L/mo time steps. Let us denote the final 1% of this signal by (Msig,t^'sig)- We then 
compute the energy Eu (^2-norm normalised by non-dimensional signal duration) of the deviation of Msig 
from its mean: 



Mo I Msig 



-(Ms: 



(34) 



where |Msig| and denote the length and mean of Msig, respectively. We choose this robust statistic instead 
of attempting to measure signal amplitude because of numerical noise in the LBM simulation. The source 
of noise in LBM is attributed to the existence of an inherently unavoidable neutral stability direction in the 
numerical scheme (see, e.g., f3l). 

We opt not to employ the "bounce-back" boundary condition used in the previous steady state study. Instead 
we will use the diffusive Maxwell boundary condition (see, e.g., |8|), which was first applied to LBM 
in I?]. The essence of the condition is that populations reaching a boundary are reflected, proportional to 
equilibrium, such that mass-balance (in the bulk) and detail-balance are achieved. The boundary condition 
coincides with "bounce-back" in each corner of the cavity. 

To illustrate, immediately following the advection of populations consider the situation of a wall, aligned 
with the lattice, moving with velocity M^aii and with outward pointing normal to the wall in the negative 
y-direction (this is the situation on the lid of the cavity with M^aii = mq). The implementation of the diffusive 
Maxwell boundary condition at a boundary site {x,y) on this wall consists of the update 



with 



fi{x,y,t + l) = 7/;*(Mwaii), = 4,7,8, 



.f2ix,y,t) + f5{x,y,t) + f(,{x,y,t) 



ft (Mwall ) + /v* (Mwall ) + /g (Mwall ) ' 

Observe that, because density is a linear factor of the equilibria, the density of the wall is inconsequential 
in the boundary condition and can therefore be taken as unity for convenience. As is usual, only those 
populations pointing in to the fluid at a boundary site are updated. Boundary sites do not undergo the 
collisional step that the bulk of the sites are subjected to. 

We prefer the diffusive boundary condition over the often preferred "bounce-back" boundary condition 
with constant lid profile. This is because we have experienced difficulty in separating the aforementioned 
numerical noise from the genuine signal at a single control point using "bounce-back". We remark that the 
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Fig. (12): Plot of energy squared, E}^ ( |34] |. as a function of Reynolds number. Re, using LBGK regularised 
with the median filter limiter with 5 — 10^^ on a 100 x 100 grid. Straight lines are lines of best fit. The 
intersection of the sloping line with the x-axis occurs close to Re = 7135. 

diffusive boundary condition does not prevent unregularised LBGK from failing at some critical Reynolds 
number 

Now, we conduct an experiment and record ( |34] | over a range of Reynolds numbers. In each case the median 
filter limiter is employed with parameter 5 = 10^^. Since the transition between steady and periodic flow 
in the hd-driven cavity is known to belong to the class of standard Hopf bifurcations we are assured that 
eI oc Re lfT2l . Fitting a line of best fit to the resulting data localises the first bifurcation in the lid-driven 
cavity flow to Re — 7135 (Fig. |(12)| i. This value is within the tolerance of Re ~ 7402 ±4% given in |28| 
for a 100 X 100 grid. We also provide a (time averaged) phase space trajectory and Fourier spectrum for 
Re = 7375 at the monitoring point (Fig. |(13)] and Fig. |(14)| i which clearly indicate that the first bifurcation 
has been observed. 

9. CONCLUSION 

In the ID shock tube we do not find any evidence that maintaining the proper balance of entropy (imple- 
menting ELBM) regularizes spurious oscillations in the LBM. We note that entropy production controlled 
by a and viscosity controlled by j3 are composite in the collision integral (|6]). A weak lower approximation 
to a would lead effectively to addition of dissipation at the mostly far from equilibrium sites and therefore 
would locally increase viscosity. Therefore the choice of the method to implement the entropic involution is 
crucial in this scheme. Any method which is not sufficiently accurate could give a misleading result. 
In the 2D lid driven cavity test we observe that implementing TRTl 13 1 or MRT ll22l with certain relaxation 
rates can improve stability. The increase in stability from using TRT can be attributed to the correction of 
the numerical slip on the boundary, as well as increasing dissipation. What is the best set of parameters to 
choose for MRT is not a closed question. The parameters used in this work originally proposed by Lallemand 
and Luo |'22'| are based on a linear stability analysis. Certain choices of relaxation parameters may improve 
stability while qualitatively changing the flow, so parameter choices should be justified theoretically, or 
alternatively the results of simulations should be somehow validated. Nevertheless the parameters used in 
this work exhibit an improvement in stability over the standard BGK system. 

Modifying the relaxation rates of the different modes changes the production of dissipation of different 
components at different orders of the dynamics. The higher order dynamics of latttice Boltzmann methods 
include higher order space derivatives of the distribution functions. MRT could exhibit the very nice property 
that where these derivatives are near to zero that MRT has little effect, while where these derivatives are large 
(near shocks and oscillations) that additional dissipation could be added, regularizing the system. 
Using entropic limiters explicitly adds dissipation locally f?|. The Ehrenfest steps succeed to stabilize the 
system at Reynolds numbers where other tested methods fail, at the cost of the smoothness of the flow. We 
also implemented an entropic limiter using MRT technology. This also succeeded in stabilizing the system 
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Fig. (13): A phase trajectory for velocity components for the signal (Msig,Vsig) at the monitoring point 
(L/16, 13L/16) using LBGK regularised with the median filter limiter with 5 — 10^^ on a 100 x 100 grid 
(Re = 7375). Dots represent simulation results at various time moments and the solid line is a 100 step time 
average of the signal. 




Fig. (14): Relative amplitude spectrum for the signal Wsig at the monitoring point (L/16, 13L/16) using 
LBGK regularised with the median filter limiter with 5 = 10^-^ on a 100 x 100 grid (Re = 7375). We 
measure a dominant frequency of o = 0.525. 



to a degree, however the amount of dissipation added is less than an Ehrenfest step and hence it is less 
effective. The particular advantage of a limiter of this type over the Ehrenfest step is that it can preserve the 
correct production of dissipation on physical modes across the system. Other MRT type limiters can easily 
be invented by simply varying the relaxation parameters. 

As previously mentioned there have been more filtering operations proposed|29|. These have a similar idea 
of local (but not pointwise) filtering of lattice Boltzmann simulations. A greater variety of variables to filter 
have been examined, for example the macroscopic field can be filtered rather than the mesoscopic population 
functions. 

We can use the Enstrophy statistic to measure effective dissipation in the system. The results from the 
lid-driven cavity experiment indicate that increased total dissipation does not necessarily increase stability. 
The increase in dissipation needs to be targeted onto specific parts of the domain or specific modes of the 
dynamics to be effective. 

Using the global median filter in the lid-driven cavity we find that the expected Reynolds number of the 
first Hopf bifurcation seems to be preserved, despite the additional dissipation. This is an extremely positive 
result as it indicates that if the addition of dissipation needed to stabilize the system is added in an appropriate 
manner then qualitative features of the flow can be preserved. 

Finally we should note that the stability of lattice Boltzmann systems depends on more than one parameter 
In all these numerical tests the Reynolds number was modified by altering the rate of production of shear 
viscosity. In particular, for the lid driven cavity the Reynolds number could be varied by altering the lid 
speed, which was fixed at 0. 1 in all of these simulations. Since the different modes of the dynamics include 
varying powers of velocity, this would affect the stability of the system in a different manner to simply 
changing the shear viscosity coefficient. In such systems the relative improvements offered by these methods 
over the BGK system could be different. 

The various LBMs all work well in regimes where the macroscopic fields are smooth. Each method has its 
limitation as they attempt to simulate flows with, for example, shocks and turbulence, and it is clear that 
something needs to be done in order to simulate beyond these limitations. We have examined various add- 
ons to well-known implementations of the LBM, and explored their efficiency. We demonstrate that add-ons 
based on the gentle modification of dissipation can significantly expand the stable boundary of operation of 
the LBM. 
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